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ABSTRACT 

The pressure amplitude distribution along the interface between a 
wedge-shaped fluid layer overlying fast and slow fluid bottoms were 
investigated using computer models based on the method of images. 
The pressure amplitude along the bottom rapidly falls off for distances 
closer to the apex than that at which adiabatic-mode theory predicts 
cut-off for the lowest mode. This distance is termed the “dump 
distance” for an acoustically fast-bottom condition. In the case ofa 
slow bottom, because the speed of sound ratio is greater than 1, the 


dump distance does not exist. To facilitate scaling the acoustic fields, 


a “scaling distance” is introduced: XO = Zino, tand where 


Cc 
He = cos! = ko = _ 8 is the angle of the wedge, c, = speed of sound 


in water, c2 = speed of sound in bottom, and A = wavelength in the 
wedge medium. Xo appears to have physical meaning as a useful 


scaling distance. 
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I. INTROD ION 


Acoustic propagation in an ocean with a sloping bottom is of 
growing concern both for theoretical reasons and because of the 
importance of knowing the performance of acoustic sensors located 
over the continental slope. This problem has been studied by many 
scientists [Ref. 1-8], but many have examined only cases of a fluid 
wedge for which both surfaces are acoustically impenetrable. 

In the case of a fluid wedge for which both surfaces are acousti- 
cally impenetrable, mathematical examination of the problem has been 
performed by Bradley and Hudimic [Ref. 1] using both image theory 
and normal-mode theory. 

For practical interest for underwater acoustics, the top boundary 
may still be approximated as an impenetrable, pressure-release sur- 
face, but the bottom surface can no longer be considered to be per- 
fectly reflecting. From a mathematic point of view, this generalization 
of the boundary conditions has an important repercussion: The wave 
equation is no longer separable and standard normal-mode theory can 
no longer be applied. 

This research will focus on four tasks: 


1. Use the models to reevaluate the fast-bottom case which was 
done by Lesesne [Ref. 9]. 


2. Provide output in graphic form for ease of interpretation. 


3. Change from the fast bottom model to a slow bottom model and 
evaluate whether a scaling factor is useful or not. 


4. Compare results to others’ research. 


eZ 


I. DEVELOPMENT 


A. DESCRIPTION OF THE THEORY 
1. Image Theory 

The approach used in this research is to determine the 
amplitude and phase of the pressure in the fluid by the method of 
images [Refs. 8-10]. In 1978, with the aid of a computer model, 
Coppens, Sanders, lonnou, and Kawamura [Ref. 11] predicted and 
measured the pressure amplitude and phase in the up-slope direction 
along the bottom of a wedge-shaped fluid layer overlying a fast fluid 
bottom. The geometry is shown in Figure l. 

Baek [Ref. 12] used the same model in 1984 to predict pres- 
sure amplitude and phase everywhere within the wedge in the up- 
slope direction. In the same year, Lesesne [Ref. 9] implemented a 
model developed by Coppens and Sanders which is not limited to up- 
or down-slope prediction. 

For the wedge-shaped duct, the images which lie on a circle 
(see Figure 2) are eventually going to appear in the region where the 
source is located and which, by the established boundary conditions, 
can contain only the source. However, following Sommerfeld [Ref. 13], 
the concept of an extended Rieman surface is introduced. The exten- 
sion of the 6 coordinate by introduction of Rieman sheets allows for an 
infinite (if need be) number of images. 

The source and each of the images radiate a spherical wave of 


the appropriate phase. The phase coherent summation of these waves 
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yields the total pressure and phase to be found at any field point in the 
wedge. 

The total pressure of all the images and the source along a 
line of constant x from the apex can be determined by the method of 


images (see Figure 2). 


_ N 
Py (x) = y (-1)INTIN/2) [g,_9 exp[jk1x cos(®n - D)]] + 
n=1 


Sn exp{jk1x cos(8n + D)] (1) 
Where 
N = INT(180/8), 
8 = the wedge angle, 
D = the receiver angle, 
G = the source angle. 
6n = (n - 1)8 + G, l= eo onier 
Oyen) oe re n=2 46.875 
Zn = R(@m) m odd 
In=3. oO. 
n= R(8m) m even 
m=2,4,6... 
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The gn-2 and gn are path parameters from the nth image and R(6,) is 
the reflection coefficient of the bottom with grazing angle 6,. (For 
details, see Ref. 12.) 

The successive reflection of sound from the surface and bot- 
tom is accounted for by multiplying the sound field radiated by each 
image by the plane wave reflection coefficients corresponding to reflec- 
tions encountered by the wave as it propagates toward the field point. 

For the fast bottom, all distances are normalized to the 
distance Xo, measured from the apex, at which the lowest mode would 
attain cut-off in the adiabatic approximation. This distance Xo is 


called the dump distance [Ref. 12] (see Figure 2). 


eS 
anB ~ 4 * sin@, * tang 
XK = wavelength in wedge medium 
8. = cos7! = = critical angle 
= wedge angle 
h = water depth at dump distance 


For this model, we use wedge angles of the form x/n where n 
is an integer and consequently no diffraction term exists [Refs. 14-15]. 
If n is not an integer, then there will be a diffracted wave which will 
always arrive at a later time than the latest image contribution 


[Ref. 14]. 
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2. Normal Mode With Ray Bundles 

V. K. Kuzentsov's [Ref. 7] normal-mode-with-ray-bundles 
approach treats the entire domain of the wedge as two zones with 
respect to the separation at Vy. where Vp, is the exit line of the nor- 
mal mode from the wedge into the half-space. Its position corre- 
sponds to the critical wedge thickness for the given normal mode 
with the non-wave guide zone r < rm. In the wave-guide zone, the 
normal modes behave as normal modes in an ideal wedge. 


3. Ray Tracing 
Tien et al. [Ref. 16] developed an analysis based on a ray 
approach. At the depth where the mode would be cut off, it is 
assumed that a ray reflects from the bottom of the wedge at the criti- 


cal angle: 


This ray is allowed to reflect from the surface and it reattains the 
bottom at a grazing angle which is the critical increased by twice the 
wedge angle. The intervening region of the wedge is filled with a 
large number of downward-traveling rays whose angles of incidence 
on the bottom are uniformly increasing between 8 and 0 + 28. Then, 
for each ray the transmission into the bottom is calculated. This 
process is repeated progressively toward the apex of the wedge. 
There is no attempt to preserve phase information, and the power 
of the beam in the bottom at each angle of depression is simply the 
sum of the powers transmitted by the appropriate rays. Calculations 
performed with the help of a computer again provide good agree- 
ment with their experimental observations as far as the angle of 
depression and beam width of the beam in the bottom were con- 
cerned, but this approach predicts that there should be no trans- 
mission of sound into the bottom until the cutoff depth is reached. 
[Ref. 17] 


In the ray-tracing method, there is no attempt to preserve 
phase information, and the power of the beam in the bottom is simply 


the sum of the power transmitted by the appropriate rays. 
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4. Par lic Equation 


Jensen and Kuperman [Ref. 4] used the parabolic equation 
approximation to study up-slope propagation. Their numerical calcu- 
lations suggested that there is negligible coupling between normal 
modes. They also obtained some very interesting graphs which show 
that the energy of a mode passing through cut-off was radiated as a 
beam into the bottom. 

The only significant limitation of this method appears to be 
that inherent in the approximations necessary for applicability of the 
parabolic equation itself. 

5. Buckingham’s Approach 

In an A. B. Wood Memorial Lecture, Buckingham [Ref. 2] dis- 
cussed acoustic propagation in a wedge-shaped ocean. Basically, he 
follows his previous papers [Ref. 18], modified to include a penetrable 
bottom. In the original paper [Ref. 18], both forward scattering and 
backward scattering exists in the acoustic field except when the 
receiver is located in the vicinity of the source. However, in the 
penetrable wedge there can be no significant backward-reflected 
component in the acoustic field [Ref. 2]. 

Buckingham also used the Rayleigh law of reflection for fast 
bottom and derived an effective wedge with a pressure-release bottom 


(see Figure 31): The Rayleigh law of reflection is 


vy . & sino - j(cos2a - cos2a,)1/2 41) 
~ g sin a + j(cos2a — cos2a,) 172 
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where g = - a is the angle from the normal, and a, is the critical 


angle. Here, we are concerned only with grazing angles less than the 
critical, where total internal reflection occurs. When a > de, there is a 
phase change on reflection, 2E, which from equation (1) can be 
expressed as 


(cos2a - cos2a¢) a 


= eal 
2 oe? an) sina 


(2) 


Buckingham then points out that similar phase change would occur if 
the ray had been reflected from a pressure-release surface a distance 


A below the actual interface (Figure 31): 


1 
~-E 
2 
(EE a 


where k = wave number in the first fluid. At the critical angle, it is 


clear from equation (2) that E = O, in which case, 


rt 
AO = dksinac 


B. DEVELOPMENT OF THE WEDGE MODEL 

Originally, the fast-bottom model was modified from a computer 
model by Coppens, which was programmed in FORTRAN by Baek and 
designed to run on the IBM 3033. The slow-bottom model resulted 
from a modification of this program. For the slow bottom, the cut-off 


condition does not exist. However, using the following formula, 
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] 
99 = cos~! a c) = velocity sound in water 


co = velocity sound in bottom 


a scaling factor is introduced which uses the same idea that was in the 


fast-bottom model but with the speed of sound ratio inverted. 


C. DEVELOPMENT OF THE GRAPHIC PROGRAMS 

At this stage, it became necessary to consider alternate forms of 
data presentation. Because of the large volume of computed values, 
numerical output had become ineffective for recognizing details. The 
graphic package Disspla [Ref. 20], available on the IBM 3033 com- 
puter, was used in conjunction with these programs. 

Graphical plots of the data are essential for a clear understanding 
of the problem. A purpose of this research was to utilize graphics to 
better determine the physical concepts necessary to interpret the 
aspects of sound propagation in a wedge-shaped ocean with a pene- 
trable bottom. 

Disspla is an option which uses “function calls” from a FORTRAN 
program to generate graphs. Disspla’s use of FORTRAN allows func- 
tion calls to be written into a program or to be written as separate 
programs. The data used to compute the graphs can, therefore, be 
generated within the FORTRAN program or read in from a data file. 
The latter method was employed since it offered more versatility. 

The size of the data file array or the number of points used did not 
prove to be a limitation on the graphs that were generated. By holding 


all but two parameters constant, a plot of pressure amplitude versus 


s, 


the remaining variable was obtained. The position of the receiver on 
the bottom was usually varied to map the pressure field for a given 
source location. 

The Disspla package supports this, but several problems had to be 
overcome. The Pocket Guide of Disspla describing this option was 
incomplete and very confusing. After comparing the various options 
available, it was decided to use the mode of graphics in which a matrix 
is used to define the data. 

At the beginning, it was necessary to determine the optimal num- 
ber of matrix points which would provide smooth and reasonable plots. 
After running through hundreds of plots, it was decided to use 100 x 
100 matrix points. 

The other important limitation was that the array containing the 
data points must be initially dimensionalized to the exact size of the 
data file array. This information was normally passed from the data file 
during the programs. This made it necessary to initialize the size of 
the data file array before the graphics program was run. 


1. Pressure-Amplitude Graphics 


a The pressure amplitude was plotted as a function of the two 
receiver coordinates. This produced a 3-D plot. The clearest 
presentation was a series of 3-D graphs stacked one after another 
(see Figures 3-A and 4-A), which contain fast- and slow-bottom 
complementary cases. (In this context, “complementary case” 
means the two cases obtained by exchanging the relative values of 
C; and co.) 


b. When the receiver is close to the source, the spatial variation of 
the pressure amplitude is smaller than the increment, producing 
artificial structure to the graph. To allow a truer view of this 
portion of the pressure amplitude, smaller portions of the field 
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can be viewed with higher resolution without increasing com- 
puter time (see Figures 5 and 6). 


c. When a fuller view was necessary, the models were modified to 
ignore the field close to the source (see Figure 7). 


2. Pressure ntour 
The second graphic program used was a contour plot. The 
following should be mentioned here: 
a The poly 3 interpolation method was used to produce the plots. 


b. The number of matrix points in this program is the same as as for 
the pressure-amplitude graphics. 


c. The increment in the pressure amplitude between contour lines 
was controlled either automatically or by the user. (In Figures 
3-B and 15-B, one can see the difference between the two 
methods.) 

In an attempt to analyze the pressure distribution near the 
source, a polar graphic program was developed. Here, the receiver 
was positioned on circles centered on the source and moved out in 
radial increments. Although this polar form satisfactorily dispersed 


the pressure field near the source, it created very serious geometrical 


distortion (see Figure 8). As a result, it was not very useful. 


21 


III. RESULTS FROM THE PLOTS 


A. DATA PRESENTATION 
1. The 3-D Plots 
Generally, the source is fixed and the pressure amplitude is 
calculated as a function of receiver position. So, for the three axes of 
the graphs: 


1. The Ro axis is the normalized receiver distance from the 
shoreline. 


2. The Yo axis is the normalized distance along the shoreline 
between the receiver and the source. 


3. The vertical axis is the pressure amplitude. 
The resolution of the plots is set by specifying a 100 x 100 


matrix of pressure calculation, which means: 


The range of Ro (or Yo) 


100 = The increment of Ro (or Yo) 


The smaller the range of Ro (or Yo), the smaller the increment of Ro 
(or Yo), the higher the resolution of the plots, and the smoother the 
plots. 
2. The Contour Plots 

The axes of the contour plot are X and Y, in which X is equal 
to Ro and Y is equal to Yo. The label values on the contour lines are the 
pressure amplitudes. The increments of the pressure amplitude and 
the resolution are also controlled by the value of CONMAK [Ref. 20] 
(see Figure 9-B). 
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There are many “tear drop” contours, as shown in Figure 4-B. 
The contours are caused by the limited resolution and disappear at 


higher resolution (see Figure 9-B). 


B. GENERAL PROPERTIES OF THE PLOTS WHICH ARE THE 
SAME FOR BOTH MODELS 


When R, is equal to an odd number, the pressure peak is directly 
under the source, as shown in Figures 10 and 13. When R] is equal to 
an even number, the pressure is a minimum under the source (see 
Figures 11 and 14). The explanation for this can be found by using 
only the source, the image character, and the first surface reflected 
path (see Figure 12). 

The time dependence of the signal at a distance d, from the 
source is 

Pi = = e-jl@t-kd}) 
and the time dependence of the signal arriving at R from the image 
characterizing the first surface reflected path is 


-A 


* 
Py = do e-j(wt-kd9) 


where dg is the distance from the image point to receiver. So, 


* A A 
= = e- t | — e-jkd; ~ — e-jkd 
By = P) + P, e-jo Res 1 do © 2 


= e-jot f-4 Jk(da-a) 
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Thus, the phase difference between the signals arriving at R from S 


and S' is 


KG Gear 


When constructive interference occurs, 


AOo= 2m, n= OFZ 


and when destructive interference occurs, 


Ao = (2a Yn on = OZ] 


The set of points for which dg - d; has a constant value defines a locus 
along which Ad is constant. The locus d2 - dj) = constant is a hyper- 
bola of revolution about the vertical axis connecting the source and the 
first surface-reflected image. The intersection of this surface and the 
plane representing the interface between the two media yields ellipti- 
cal curves close to the source, which evolve into hyperbolic curves far 
from the source. 

These features are observed in the pressure contour plots for the 
various cases. 

The angle of intromission for which the reflection coefficient is 
zero will exist in the fast-bottom case for p; > pg or the slow-bottom 
case for 91 <p2. Figures 17 and 18 are for cases possessing angles of 
intromission; these appear to exhibit no significant pressure changes 


compared to the non-intromission conditions. 
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C. DIFFERENCE BETWEEN BOTH MODELS 

With the fast bottom, the pressure pattern is much more irregular 
and the pressure decays with distance from the source less rapidly 
because, with the fast bottom, more images are totally reflected from 


the bottom. 
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IV. EVALUATION 


The effects of changing the properties of the wedge and the 
source position on the prediction of the pressure field were evaluated. 


As a base-line case, the following parameters were used: 


B a leOn 
= 5° 
ie = Q° 
R} = 40 
= = 09 
c)}/c9 = 0.9 (fast bottom) 
C1, co = 1/0.9 (slow bottom) 
a/ko = 26.0001 


Three-dimensional depictions of the pressure fields for these 
conditions as functions of Rg and Yo are shown in Figures 3 and 4 for 
both fast- and slow-bottom cases. The corresponding contour plots 
are included. 

Both models were tested by varying one parameter and holding 
the others constant. Comparison of pressure-contour plots as well as 
three-dimensional plots was then used to determine the effect of each 


parameter. 


A. SPEED-OF-SOUND RATIO, DENSITY RATIO, AND ATTENUATION 
The base-line case was first examined by changing the speed-of- 


sound ratio. The positions and number of crests exhibited simple and 
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clear patterns for the slow-bottom model (see Figure 19) and more 
complex patterns for the fast-bottom model (see Figure 20) as this 
ratio was changed from 0.9 (or 1/0.9) to 0.5 (1/0.5). The density ratio 
was then changed down to a value of 0.5. The positions and number of 
crests show only minor variation (see Figure 21). This suggests that 
the major features of the sound field in the wedge are relatively insen- 
sitive to minor variations in the characteristics of the bottom. 

Finally, the value of a/kg was increased from 0.0001 to 0.1. The 
pressure-pattern did not change significantly but, as expected, the 
amplitudes of the crests showed sharp decreases as the attenuation 


increases. 


B. RECEIVER ANGLE 

The reader should be aware of the fact that the program con- 
tained in the appendix was modified as simply as possible, to insure 
that pressure was calculated infinitesimally above the bottom. If it is 
desired to calculate the pressure for receiver angles not essentially 


zero, the leading statements and do loops must be modified. 


C. WEDGE ANGLE 

As the wedge angle decreases, the dump distance increases, the 
number of images increases, and therefore the CPU time increases. 
For example, for B = 10°, the CPU times for fast-bottom and slow -bot- 
tom cases are about 520 seconds and 480 seconds, respectively. If B = 
7°, the CPU time increases to about 900 seconds (See Figures 22 and 


23). 
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D. SOURCE ANGLE 

Changing the distance of the source from the apex (but still with 
G = 8/2) had a great effect on the position of the crests and troughs. 
However, for the source at large distances from the apex. scaling the 
receiver position restores the original patterns. For example, if the 
source is moved outward to twice its original position, a plot with the 
length scales divided by two produces the same graphic presentation 
(see Figures 3 and 24 or 4 and 25). This is consistent with Lesesne’s 


observation. 


E. LLOYD’S MIRROR 

Figures 29 and 30 show the field due to the source and the first 
image point only. The source and the first image are of equal strength. 
The first image is coherent with the source but nm out of phase (see 
Figure 12). That produces the Lloyd’s mirror. The pressure 
distribution observed is identical for both fast- and slow-bottom 
programs when they were modified to have the same apparent speed 
of sound in the wedge. These figures show that major features of the 
pressure field are due mainly to the source and first image when the 


receiver is close to the source. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


The theory of 3-D acoustic propagation in a penetrable wedge 
presented in this research is based on the method of images which 
includes the fast- and the slow-bottom models. One of the research 
goals was to test whether the scaling factor method can apply to the 


slow-bottom case. 


A. CONCLUSIONS 

For both the fast bottom and the slow bottom, when the receiver 
is near the source, strong spatial oscillations in pressure are present. 
These variations are mainly due to interference between the source 
and the first surface-reflected image. 

The primary result of the fast-bottom model is the importance of 
interference evident in the form of strong spatial oscillation. When 
the bottom mismatch increases, more interference occurs, because 
more images contribute to the field. 

Generally, investigation of this model showed that the scaling fac- 
tor worked. 

One observation is in violation of physical intuition—when the 
bottom mismatch increases, more interference should occur. Figures 
26, 27, and 28 show the opposite, with less interference as the bot- 
tom mismatch increases. 

The results for the fast-bottom model are more complex than for 


the slow-bottom case because the cut-off condition only exists in the 
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fast-bottom model, which yields more images which correspond to 
perfect reflections for the bottom. 

The method of images does not yield an easy decomposition of the 
sound field into a form that can be closely associated with normal 
modes. On the other hand, the method of images provides intrinsi- 
cally the phase interference effects associated with the superposition 
of effects that would be more difficult to extract from the adiabatic 
normal mode. The intrinsic inclusion of absorption in the substrate is 
an attractive and realistic feature of the image approach. In addition, 
the image approach appears to avoid the generation of caustics at 
turning points which are inextricably linked with the results of 


unmodified ray tracing approaches. 


B. RECOMMENDATIONS FOR FURTHER INVESTIGATION 
1. For the case of a real fast-bottom, identification of the normal 
modes of the system and investigation of the possibility of normal 
mode coupling is worth study. 


2. Using the batch system to investigate the small wedge angle is 
worth study. 


3. Changing the position and depth of the receiver is worth study. 


4. For the slow bottom, experiments should be undertaken to com- 
pare with the theoretical predictions. 
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APPENDIX A 
COMPUTER PROGRAMS 


FILE: 3LS4 FORTRAN A1 


DEBUG UNIT(6),SUBCHK 
END DEBUG 


C CROSS SLOP TRANSMISSION IN THE WEDGE WITH PENETRABLE LOSSY SLOW 
C BOTTOM (NO GEOMETRICAL APPROXIMATIONS ) 


@ DIMENSIONZ(M,M) (21,11), (Y0,R2)--(L,A) 

DIMENSION Z(101,101) 

COMMON WORK(50000) 

INTEGER A\lI,JJ,N1,$1,S$2,1E,N, 

REAL*4 RR1,R2,Y0,Y1,YOSTA,YOINC,R2INC,R2STA,RZ,RZ1,M.J, 
1 B,C(30),CC,C2,D,01,D2,E(30),E1,E2,F(30),F1,-2,G,H,PI,FI, 
1 P1,P2,P(125,130),PH,Q1,Y2,R3,R8(80),R9(80), 
1 §(30),T,T4,T6,T1(80), TQQ,TQQ1,TQAQ2,TQAQ3,W,W0,W1,XL,Y, 
1 21,Z2,Z3,24,25,Z6,Z,EL,DIST 

CALL SETIME 

P| = ARCOS(-1.0) 

INPUT PARAMETERS 

N1 = # OF IMAGINE POINTS 

B = WEDGE ANGLE (DEG) 

G = SOURCE ANGLE (DEG) 

D = RECEIVER ANGLE (DEG) 

R1 = SOURCE DISTANCE (IN DUMP DISTANCES) 

R2 = RECEIVER DISTANCE (IN DUMP DISTANCES) 

YO = SHORE DISPLACEMENT (IN DUMP DISTANCES) 

D1 = RHO1/RHO2 

SC SOC 

XL = ALPHA/K> 

A = # OF RECEIVER POSITIONS IN R2 DIRECTION 

L = # OF RECEIVER POSITIONS IN YO DIRECTION 

R2INC = INCREMENTAL INCREASE FOR R2 

YOINC = INCREMENTAL INCREASE FOR YO 


RERAREKRAEEEAEAEERAEARAEAKAETAEAERAEAAEKERRAERERERERKERREREREKeEReEEREREREEReEKReKeeeekeeeeekekeeeeekekeeeaeteeaaeeteaeteae 


teteeeratsessasteseareataneasanaeraeeeaal AT A VALUES ##**#**teeteeneeenneanne AeRHERARREREREREERREEE 


CALL TEK618 
CALL COMPRS 
CALL PAGE(11,8.5) 
B=10. 

G = 1.0°B/4. 


DQ AAQAQDQIDIAADAAAAAAAAAAQAIAADYD 


31 


G=5 
O10 
Beano 
sh = 65) 
YO =0. 
B1=9 
Dae 
FOR SLOW-BOTTOM CASE 
CC = 1109 
FOR FAST-BOTTOM CASE 
Ce=00 
66-5 
XL = .0001 
A= 100 
4100 
R2INC =.1 
YOINC = .1 
Cc teeteetrerensereeeesereseeeeseree ene li A| A] PROGRAM ttt tttteteeeeneeeereeenerereeateterenteeneeene 
YOSTA = YO 
R2STA = R2 
W =0.0 
WRITE(6,300) B,G,D 
300 FORMAT (3X,’B = '.F5.2,3X,’G = ’.F5.2,3X,’D = ',F5.2,3X) 
310 FORMAT (1X,F5.2,1X,F5.2,1X,F5.2) 
WRITE(6,400) R1 
400 FORMAT (3X,'SOURCE RANGE = ',F8.2) 
410 FORMAT (1X,F6.2) 
WRITE(6,500) D1,CC,XL 
500 FORMAT (3X,’DENSITY RATIO = ’,F8.6,3X,'SPEED RATIO = ',F8.6,3X, 


OQ (O55. °@ 


: ‘ALPHA/K2 = ’,F8,6) 
510 FORMAT (1X,F6.4,1X,F6.4,1X,F6.4) 
WRITE(6,600) 
600 FORMAT (5X,’42’,6X,’Y0’,10X,'P’) 
C WRITE(22,700) R2,Y0 


700 FORMAT (1X,F6.2,1X,F6.2) 
N1 = (180/B + .00001) 
T6 = 180/PI 
B=B/T6 
G=G/T6 
C2=66.2 


HM=0 
HH=1 
D =HM*B/A 
DD =B/A 
V = 2*B/10 
IF(D.LT.V) GO TO 110 
IF(D.GE.V) GO TO 120 
110 DX =D 
C THIS DO LOOP CALCULATES THE TMETA(N) AND THE IMINAGE SLANT 
C RANGES R8(N) AND R9Q(N). 
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DO 55 HH = 1.10 
D = DX+(HH-1)*B/(10*A) 
DD = B/(10*A) 


120 Si 


T4 = P1/(2*TAN(ARCOS(1/CC))* TAN(B)) 
TQQ = TAN(B) 
FOR SLOW-BOTTOM CASE 
TQQ, = ARCOS(1/CC) 
FOR FAST—BOTTOM CASE 
TQQ, = ARCOS(CC) 
TQQ2 = TAN(TQQ1) 
TQQ3 = 2.*TQQ2*TQQ 
T4 = PI/TQQ3 
Ol = SOR 2.0} 
LP1X=L+1 
IAP 1X=A+1 
DO 20 M=1,LP1X 
DO 10 J=1,IAP1X 
D2 = (Y0**2)+(R1**2)+(R2**2) 
R3 = 2*R1*R2 
DO 30 N=1,N1 
IF(S1.GT.0) T1(N)=FLOAT(N-1)*B+G 
IF(S1.LT.0) T1(N)=FLOAT(N)*B-G 
$1i=-S$1 
R8(N) = SQRT(D2-R3*COS(T1(N)-D)) 
R9(N) = SQRT(D2-R3*COS(T1(N)+D)) 
CONTINUE 


Piolo 
P2 = 0.0 
DO 40 N=1,N1 
S2 = (-1)**(IFIX(N/2+.0001)) 
W1 = 2*C2*XL 


oS 88 56 a Se ee oe es ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee 


11 = IFIX((N-.9999)/2) 
DO 41 1=1,11 
S(I) = ABS(R1*SIN(T1(N)-2*FLOAT(I)*B) 
+R2*SIN(2*FLOAT(I)*B-D))/R8(N) 
C(I) =SQRT(1-(S(I)**2)) 
T = S(I)/D1 
WO = (-C2+(C(I)**2)) 
Y = SQRT((W0**2)+(W1*+2)) 
Z(M,J) = ABS(WO) 
IF(Y,LE,Z(M,J)) Y = Z(M,J) 
Y1 = Q1*SQRT(Y+WO) 
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Y2 = -Q1*SQRT(Y-WO) 
E() = (T**2)-(¥2**2)-((¥1**2))/(((T-¥2)**2)4(Y 19#2) 
F(I) = (2*Y¥1*T)/(((T-¥2)**2)4(¥1*42)) 


Aq CONTINUE 
CC anne ntact nn ete ca 
z ~---eeeeennee PRODUCT OF REFL COEFFS ALONG NTH UPPER PATH--------------------- 
E1=1.0 
Fi=0.0 
IF(N.LE.2.0) GOTO 422 
DO 42 1=1,11 
E2 = E1°E(|)-F1*F(I) 
F2 = F1*E(I)-E1*F(I) 
E1=E2 
F1=F2 
42 CONTINUE 
ce AAAEEEEEEEEEEEEEEEEEEEE EERE EERE EREREEREREEEEE EERE RE EEE Eee 
Cc FOR SLOW-BOTTOM CASE 
422 T = Ta*Ra(N) 
c FOR FAST-BOTTOM CASE 
422 T = T4*Re(NV/CC 
G SERA AEREEEEEREEEEEEEEEEEEE HEE EEEEEEEEEEEEEEEEEEEEEEEE REE EE Eee 


P1 = P1 + S2*(E1*COS(T)+F1*SIN(T))/48(N) 
P2 = P2 + S2*(F1*COS(T)+E1*SIN(T))/48(N) 


S(l) = ABS(R1*SIN(T1(N)-2*FLOAT(I-1)*8) 
+ R2*SIN(2*FLOAT(I-1)*B+D))/R9(N) 
C(I) = SQRT(1-(S(1)**2)) 
T =S(IV/D1 
WO = (-C2+(C(1)*#2)) 
Y = SQRT((WO0**2)+(W1**2)) 
Z(M,J) = ABS(WO) 
IF(Y.LE.Z(M,J)) Y = Z(M,J) 
Y1 = Q1*SQRT(Y+WO) 
Y2 = -Q1*SQRT(Y-WO) 
E(l) = ((T**2)-(¥2**2)-(Y 1**2))/(( (Ty vine 2 ental ie a) 
F(I) = (2*Y¥1*T)/(((T-¥2)**2)+(Y 1#*2)) 
c AN2 = TAN(F(I)/E(1))*T6) 
43 CONTINUE 


F1=0. 

DO 44 1=1,11 
E2 = E1*E(I)-F1*F(I) 
F2 = F1*E(1)+E1*F(I) 
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Ee? 
F1=F2 
44 CONTINUE 
G rrrrrrrrrrr et ee tir errr rrrrrrrrrrrrerrrrcercrrcrcrcercrcercrercerreer es | 
€ FOR SLOW-BOTTOM CASE 
T = T4*Ro(N) 
C FOR FAST-BOTTOM CASE 
T = T4*Ro(N)CC FAST BOTTOM 
G eee ARAARAAERERAAARARARAEEAEREAEEARAER ERA RARE ERE EERE EERE REE 
P1 = P1 + S2*(E1*COS(T)+F1*SIN(T))/R9(N) 
P2 = P2 + S2*(F1*COS(T)-E1*SIN(T))/R9(N) 
40 CONTINUE 
Z(M,J) = SQRT((P1**2)+D(P2**2))*R1 
IF(W.LT.Z(M,J)) W=Z(M.J) 
WRITE(6,200) R2,Y0,Z(M.J) 
200 FORMAT (3X,F7.3,3X,F7.2,3X,F9.5,3X,F9.5) 


Se RHEE AEREREEEREREAKRERERERAARKEEKREEARKREAREREReEAHReRARARe a HeAR AT HERARARHEAE EKER E 


CG RHHKKEKREKEERERKERERERREREREEEE RAE REAR AEA AREA AHReERE RARER AAA eERE EEE REET 


R2 = R2+R2iNC 


C ReEKeRAREERRREEEREARERERERERERREREREREREREEREREREKReE EERE KREREREEREKAREKREKREREREReRARRAERERARKRERE RE RRR KERR RHA EE 

S DIST = SQRT ((J-41)**2) + (M-1)*#2) 

c IF (DIST.LE.20.) Z(M,J) = 0. 

jG HEREKAKKEREEEREEHEREEEREREKEEERAEAREEAREEEHRERTEHREREREEEKRERHAEEEEERREREKEEKEREEKRE RE HRERHRRERRERARE RARE REAR E KR AREESE 
Roe. 


20 CONTINUE 
CALL SURF3D(Z,NX,NY) 
YO =0. 

55 CONTINUE 

33. CONTINUE 


STOP 
END 
nna Es a eee oe me eee ieee 
SUBROUTINE SURF3D(Z,NX,NY) 
i a 


DIMENSION Z(101,101) 
DIMENSION R2Z(101,L101) 


COMMON WORK(50000) 
are 
C CHECK DATA STATEMENTS TO ET AXIS LIMITS: 
Sk 0d a Ss oo eee So eee 

DO 1 1-1,101 

DO 1 J-1,101 

1 R2Z(J,l)=Z(I,J) 

iO ew le ns So a a ie a ee 
CO nenne mene nnn een PRODUCE CONTOUR PLOT--—--—--—--=--=--=--= 7+ 
SO eee ae ene 


CALL PAGE (11,8.5) 
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CALL SHDCHR(.90,1,,003,1.) 
CALL AREA2D(9.5,6) 
Cc CALL NOBRDR 
CALL BLOWUP(.85) 
CALL XNAME ('XXX’,1) 
CALL YNAME(‘YYY’,1) 
CALL HEADIN(‘CONTOUR PLOT$’,100,1.,1) 
CALL GRAF(5.,1.,15.,0.,1., 10.) 
CALL FRAME 
CALL BCOMON(50000) 
CALL CONMAK(RZ,101,101,.5) 
CALL CONLIN(0,'SOLID’,'LABELS’,1,1) 
CALL CONANG(90.) 
CALL POLY3 
CALL CONTUR(2,’LABELS’,DRAW’) 
CALL ENDPL(0) 


2 
CG, PRODUCE SURFACE PLOT-—-—-- 
6 NID en ce MENS Tas i DE 
CALL PAGE(11,8.5) 
CALL AREA2D(9.5,6) 
CALL MX1ALF(’STANDARD’,’*’) 
CALL MX2ALF(’L/CSTD’,'+’) 
CO Sc cee eee 


C LABEL AXES: 
CALL X3NAME(’R2 A+XIS*$',100) 
CALL Y3NAME(‘Y0 A+X!S*$',100) 
CALL Z3NAME(‘PRESSURE A+XIS*$’,100) 


Ce ates eesti ola peo te 
CALL VOLM3D(5.,5.,5.) 
CALL YINTAX 
CALL XINTAX 
CALL ZINTAX 
CALL ZAXANG(0.) 
CALL YAXANG(0.) 
CALL XTICKS(5) 
CALL YTICKS(5) 
CALL HEADIN (‘SLOW BOTTOM’,11,1.5,1) 
CALL HEADIN (‘FAST BOTTOM’, 11,1.5,1) 
CALL MESSAG (’C.P.U. = SECONDS$',100,6.5,0.) 
CALL MESSAG (‘VIEPTS IS 250,350,90$',100.6,6.5) 
CALL MESSAG (‘SURMAT IS RZ,1,101,101,101$',100,6,6) 
CALL MESSAG (‘SURVIS IS TOP$’,10,6,5.5) 
CALL MESSAG ('B(W.A.)-10.,G(S.A)-5.,D(R.A)=0 $',100,6,5) 
CALL MESSAG (‘RHO1/RHO2 =.9,ALPHA/K2=.0001 $',100,6.4.5) 
CALL MESSAG ('C1/C2 -1/.9,R1=10.,42=5.,Y0=0.$',100,6,4) 
CALL MESSAG (‘R2INC =.1,YOINC =.1,CMK =.5$',100,6,3.5) 
Caen tee ah Scns 


CALL VUABS(250.,350.,90.) 
CALL SURVIS(‘TOP’) 
CALL GRAF2D (5 1p i520 ceo) cmon 
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ALL BLSUR 
c CALL SURMAT(RZ,XX,A+1,XX,L+1) 
CALL SURMAT(RZ,1,101,101,101,0) 
ED Pe ee ee a oo asa oon a actin amen oni ans ee ase nsneeen ce ewe ee eenee——--- 
CALL ENDPL(0) 
CALL DONEPL 
CALL GETIME(IET) 
EL = IET « .000026 
WRITE (4,2) EL 
2 FORMAT(20X, TIME = ',E16.6,,SECONDS.') 
RETURN 
END 
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Figure 2 
The Geometry of Image Solution 
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Figure 4-B 


Equal-Pressure Contours on the Bottom of the Wedge 


for the Indicated Values of Parameters 
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Figure 8-B 


Equal-Pressure Contours on the Bottom of the Wedge 


for the Indicated Values of Parameters 
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Figure 12 


Interference Pattern from a Point Source and an Image Point 
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Pressure on the Bottom of the Wedge as Calculated by the Image Model 
for the Indicated Values of Parameters 
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Figure 22-B 


Equal-Pressure Contours on the Bottom of the Wedge 


for the Indicated Values of Parameters 
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Figure 30-B 


Equal-Pressure Contours on the Bottom of the Wedge 
for the Indicated Values of Parameters 
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